################################################
# Analysis for:
#
# Warner, Zach, Garrett N. Vande Kamp, and Soren Jordan.
# ``Conditional Relationships in Dynamic Models.'' 
#  Political Science Research & Methods
#
# Required files: WVKJ-MC-simulations.csv
#                 main-cavari-psq.dta
#                 interp_urdf.R (for interpretation)
#
################################################

# Set your working directory here

library(dplyr)
library(ggplot2)
library(grid)
library(haven)
library(lmtest)
library(dynamac)
library(urca)
library(tseries)
library(polywog)
library(MASS)

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  plots <- c(list(...), plotlist)
  numPlots = length(plots)
  if (is.null(layout)) {
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                     ncol = cols, nrow = ceiling(numPlots/cols))
  }
  if (numPlots==1) {
    print(plots[[1]])
  } else {
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
    for (i in 1:numPlots) {
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

# ur.df output is annoying. interp_urdf.R helper function for interpretation
#  Author: Hank Roark. 
#  https://gist.github.com/hankroark/968fc28b767f1e43b5a33b151b771bf9
source("interp_urdf.R")

# These parameters are calso called in the sim analysis below the data creation.
# everything stationary with rho = 0.5
rho.x <- c(0.1, 0.5, 0.9) 
rho.z <- c(0.1, 0.5, 0.9) 
rho.y <- c(0.1, 0.5, 0.9) 
nobs <- c(16, 31, 51, 101, 501)  # add 1 because we'll drop the first value due to lagged values being NA

#set up the master dataframes for conditions to vary
exper <- expand.grid(nobs = nobs,
                     rho.x = rho.x,
                     rho.z = rho.z,
                     rho.y = rho.y)

# true values not varying
alpha.0 <- 0   # intercept
alpha.1 <- .5  # l.y
beta.0 <- -.3   # x
beta.1 <- .3   # l.x
beta.2 <- .2  # z
beta.3 <- -.7  # l.z
beta.4 <- .5  # x_z
beta.5 <- -.25  # l.x_z
beta.6 <- .1  # x_l.z
beta.7 <- -.4  # l.x_l.z


df <- read.csv("WVKJ-MC-simulations.csv")

df$var <- as.character(df$var)
df$model <- as.character(df$model)
df$dgp <- as.character(df$dgp)


### Reset the seed so we get the same results whether or not you've run the MC sims
set.seed(983405201)

# Add pretty plot labels
df$t <- df$nobs
df$t[df$t == 16] <- "T = 15"
df$t[df$t == 31] <- "T = 30"
df$t[df$t == 51] <- "T = 50"
df$t[df$t == 101] <- "T = 100"
df$t[df$t == 501] <- "T = 500"

df$t <- factor(df$t, levels = c("T = 15", "T = 30", "T = 50", "T = 100", "T = 500"))
table(df$nobs, df$t)

df$rho.y.fac <- recode_factor(df$rho.y, `0.1` = "rho[y] == 0.1", 
															`0.5` = "rho[y] == 0.5",
															`0.9` = "rho[y] == 0.9")

df$rho.x.fac <- recode_factor(df$rho.x, `0.1` = "rho[x] == 0.1", 
															`0.5` = "rho[x] == 0.5",
															`0.9` = "rho[x] == 0.9")

df$rho.z.fac <- recode_factor(df$rho.z, `0.1` = "rho[z] == 0.1", 
															`0.5` = "rho[z] == 0.5",
															`0.9` = "rho[z] == 0.9")

df$model.fac <- df$model
df$model.fac[df$model.fac == "full"] <- "Full interaction"
df$model.fac[df$model.fac == "restricted"] <- "Restricted"
df$model.fac[df$model.fac == "noint"] <- "No interaction"
df$model.fac <- factor(df$model.fac, levels = c("No interaction", "Restricted", "Full interaction"))


# Finally, define a df with reasonable ts
df.ts <- df[which(df$nobs %in% c(51, 101, 501)), ]








##############################################################################################################################
# Problem #1 -- bias in ECR. We should see (1) worse estimates, (2) lower coverage rates
##############################################################################################################################
# alpha estimates
a1 <- df.ts[which(df.ts$var == "alpha.1"),]
a1$true <- ifelse(a1$dgp == "spurious", a1$rho.y, alpha.1)
a1$location <- as.numeric(factor(a1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
a1$location[which(a1$model == "noint")] <- a1$location[which(a1$model == "noint")] - 0.1
a1$location[which(a1$model == "full")] <- a1$location[which(a1$model == "full")] + 0.1

# coverage
cov1 <- df.ts[which(df.ts$var == "cov.alpha.1"),]
cov1$location <- as.numeric(factor(cov1$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
cov1$location[which(cov1$model == "noint")] <- cov1$location[which(cov1$model == "noint")] - 0.1
cov1$location[which(cov1$model == "full")] <- cov1$location[which(cov1$model == "full")] + 0.1
cov1$true <- rep(.95, nrow(cov1))

# alpha plot
fig1.1 <- ggplot(a1, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = rho.y, yend = rho.y), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "ECR estimates", x = "Data-generating process",
       y = expression(paste("Mean ", hat(alpha)[1], " for each set of 500 replicates", sep = ""))) 
       
# coverage plot
fig1.2 <- ggplot(cov1, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "ECR coverage", x = "Data-generating process",
       y = expression(paste(hat(alpha)[1], " coverage probability for each set of 500 replicates", sep = "")))

######### Figure 1 #########
pdf("fig1.pdf", width = 6, height = 8)
multiplot(fig1.1, fig1.2, cols = 1)
dev.off()

######### Table 1 #########
a1b <- df.ts[which(df.ts$var == "bias.alpha.1"), ]
a1b.tab <- matrix(c(mean(a1b$median[which(a1b$model == "noint" & a1b$dgp == "spurious")]),
                mean(a1b$median[which(a1b$model == "restricted" & a1b$dgp == "spurious")]),
                mean(a1b$median[which(a1b$model == "full" & a1b$dgp == "spurious")]),
                mean(a1b$median[which(a1b$model == "noint" & a1b$dgp == "no_interaction")]),
                mean(a1b$median[which(a1b$model == "restricted" & a1b$dgp == "no_interaction")]),
                mean(a1b$median[which(a1b$model == "full" & a1b$dgp == "no_interaction")]),
                mean(a1b$median[which(a1b$model == "noint" & a1b$dgp == "restricted")]),
                mean(a1b$median[which(a1b$model == "restricted" & a1b$dgp == "restricted")]),
                mean(a1b$median[which(a1b$model == "full" & a1b$dgp == "restricted")]),
                mean(a1b$median[which(a1b$model == "noint" & a1b$dgp == "full_interaction")]),
                mean(a1b$median[which(a1b$model == "restricted" & a1b$dgp == "full_interaction")]),
                mean(a1b$median[which(a1b$model == "full" & a1b$dgp == "full_interaction")])),
                ncol = 4)
dimnames(a1b.tab) <- list(c("No interaction model", "Restricted model", "General model"),
                      c("'LDV'", "No interaction", "Restricted", "Full interaction"))

write.table(round(a1b.tab, 2),
            file = "tab1.txt",
            sep = "&",
            row.names = TRUE,
            col.names = TRUE,
            quote = FALSE)

### These are the in-paragraph interpretations
## mean coverage probabilities (these are from the in MS fig 1: t = 50, 100, 500)
summary(cov1$mean[which(cov1$model == "full")])
summary(cov1$mean[which(cov1$model == "noint")])
summary(cov1$mean[which(cov1$model == "restricted")])
# mean coverage probabilities, excluding the case it's most favorable to general model
summary(cov1$mean[which(cov1$model == "full" & cov1$dgp != "full_interaction")])
summary(cov1$mean[which(cov1$model == "restricted" & cov1$dgp != "full_interaction")])
summary(cov1$mean[which(cov1$model == "noint" & cov1$dgp != "full_interaction")])






##############################################################################################################################
# Problem #2 -- overfitting
##############################################################################################################################
# restrict to RMSE
rmse <- df.ts[which(df.ts$var == "RMSE"), ]
rmse$location <- as.numeric(factor(rmse$dgp, 
                                   levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
rmse$location[which(rmse$model == "noint")] <- 
  rmse$location[which(rmse$model == "noint")] - 0.1
rmse$location[which(rmse$model == "full")] <- 
  rmse$location[which(rmse$model == "full")] + 0.1

######### Figure 2 #########
pdf("fig2.pdf", width = 6, height = 8)
ggplot(rmse, aes(x = location, y = median)) +
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  facet_grid(vars(t)) + 
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
  scale_x_continuous(breaks = seq(1, 4, 1),labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(x = "Data-generating process", y = "RMSE") #+
 # theme(text = element_text(family = "minpro"))
dev.off()

### Global mean RMSEs
mean(rmse$mean[which(rmse$model == "full")])
mean(rmse$mean[which(rmse$model == "restricted")])
mean(rmse$mean[which(rmse$model == "noint")])

######### Table 2 #########
tab2 <- matrix(c(mean(df.ts$median[which(df.ts$var == "var.alpha.0" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.alpha.0" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.alpha.0" & df.ts$model == "full")]),
                mean(df.ts$median[which(df.ts$var == "var.alpha.1" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.alpha.1" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.alpha.1" & df.ts$model == "full")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.0" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.0" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.0" & df.ts$model == "full")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.1" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.1" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.1" & df.ts$model == "full")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.2" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.2" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.2" & df.ts$model == "full")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.3" & df.ts$model == "noint")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.3" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.3" & df.ts$model == "full")]),
                NA, # no interaction model doesn't estimate an interaction term (duh)
                mean(df.ts$median[which(df.ts$var == "var.beta.4" & df.ts$model == "restricted")]),
                mean(df.ts$median[which(df.ts$var == "var.beta.4" & df.ts$model == "full")])),
                ncol = 7)
dimnames(tab2) <- list(c("No interaction model", "Restricted model", "General model"),
                      c("alpha_0", "alpha_1", "beta_0", "beta_1", "beta_2", "beta_3", "beta_4"))

write.table(round(tab2, 2),
            file = "tab2.txt",
            sep = "&",
            row.names = TRUE,
            col.names = TRUE,
            quote = FALSE)

### size of the differences
## in-text interpretations
b4 <- df.ts[which(df.ts$var == "beta.4"), ]

mean(b4$hi[which(b4$model == "full")] - b4$lo[which(b4$model == "full")])
mean(b4$hi[which(b4$model == "restricted")] - b4$lo[which(b4$model == "restricted")])






##############################################################################################################################
# Problem #3 -- size. combine all T = 15 into one plot
##############################################################################################################################
df.15 <- df[which(df$nobs %in% c(16)), ]

a1.15 <- df.15[which(df.15$var == "alpha.1"),]
a1.15$true <- ifelse(a1.15$dgp == "spurious", a1.15$rho.y, alpha.1)
a1.15$location <- as.numeric(factor(a1.15$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
a1.15$location[which(a1.15$model == "noint")] <- a1.15$location[which(a1.15$model == "noint")] - 0.1
a1.15$location[which(a1.15$model == "full")] <- a1.15$location[which(a1.15$model == "full")] + 0.1

# coverage
cov1.15 <- df.15[which(df.15$var == "cov.alpha.1"),]
cov1.15$location <- as.numeric(factor(cov1.15$dgp, levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
# Dodge the estimated models from each other
cov1.15$location[which(cov1.15$model == "noint")] <- cov1.15$location[which(cov1.15$model == "noint")] - 0.1
cov1.15$location[which(cov1.15$model == "full")] <- cov1.15$location[which(cov1.15$model == "full")] + 0.1
cov1.15$true <- .95

########## Main MS ##########
# alpha plot
fig3.1 <- ggplot(a1.15, aes(x = location, y = median)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = rho.y, yend = rho.y), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") + 
  geom_segment(aes(x = 2.75, xend = 3.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = alpha.1, yend = alpha.1), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
   theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "ECR estimates: T = 15", x = "Data-generating process",
       y = expression(atop(paste("Mean ", hat(alpha)[1], " for", sep = ""), paste("each set of 500 replicates", sep = ""))))
       
# coverage plot
fig3.2 <- ggplot(cov1.15, aes(x = location, y = mean)) +
  theme_bw() +
  # facet_grid(vars(t)) + 
  geom_segment(aes(x = .75, xend = 1.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 1.75, xend = 2.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 2.75, xend = 3.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = 3.75, xend = 4.25, y = 0.95, yend = 0.95), color = "gray30", linetype = "dashed") +
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
  # coord_cartesian(ylim = c(0, 1)) +
  scale_x_continuous(breaks = seq(1, 4, 1), labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "ECR coverage: T = 15", x = "Data-generating process",
	   y = expression(atop(paste(hat(alpha)[1], " coverage probability for", sep = ""), paste("each set of 500 replicates", sep = ""))))
	   
# restrict to RMSE
rmse.15 <- df.15[which(df.15$var == "RMSE"), ]
rmse.15$location <- as.numeric(factor(rmse.15$dgp, 
                                   levels = c("spurious", "no_interaction", "restricted", "full_interaction")))
rmse.15$location[which(rmse.15$model == "noint")] <- 
  rmse.15$location[which(rmse.15$model == "noint")] - 0.1
rmse.15$location[which(rmse.15$model == "full")] <- 
  rmse.15$location[which(rmse.15$model == "full")] + 0.1

fig3.3 <- ggplot(rmse.15, aes(x = location, y = median)) +
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill=NA, size=1, linetype = "solid")) +
  theme_bw() +
  #facet_grid(vars(t)) + 
  geom_point(aes(shape = factor(model)), size = I(1.8), alpha = I(0.3)) +
  theme(legend.position = "none") +
 # coord_cartesian(ylim = c(.98, 1.88)) +
  scale_x_continuous(breaks = seq(1, 4, 1),labels = c("LDV", "No interaction", "Restricted", "Full interaction")) +
  labs(title = "RMSE: T = 15", x = "Data-generating process", y = "RMSE") #+
 # theme(text = element_text(family = "minpro"))

######### Figure 3 #########
pdf("fig3.pdf", width = 6, height = 8)
multiplot(fig3.1, fig3.2, fig3.3, cols = 1)
dev.off()





##############################################################################################################################
# What happens with invalid restrictions?
##############################################################################################################################
### plot estimates for beta_4
b4 <- df.ts[which(df.ts$var == "beta.4" & df.ts$dgp == "full_interaction"), ]
# Different here too: T changes the number of conditions
#  We're basically just randomly ordering them. repeat twice so that the experimental conditions are in the same order
b4$xloc <- rep(sample(seq(1:(dim(b4)[1]/2))), 2) #sample() just randomly orders the 81 experimental conditions
# clean up names for labels
b4$model[b4$model == "restricted"] <- "Restricted model"
b4$model[b4$model == "full"] <- "General model"
b4$model <- factor(b4$model, levels = c("Restricted model", "General model"))

######### Figure 4 #########
pdf("fig4.pdf", width = 6, height = 3)
ggplot(b4, aes(x = xloc, y = median)) +
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  theme_bw() +
  facet_grid(. ~ model) +
  geom_point() +
  geom_hline(aes(yintercept = .5), color = "gray30", linetype = "dashed") +
  geom_segment(aes(x = xloc, xend = xloc, y = lo, yend = hi)) +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank()) +
  labs(title = NULL, x = NULL, y = "Mean estimates with 95% quantiles")# +
#  theme(text = element_text(family = "minpro"))
dev.off()





##############################################################################################################################
# In-text replication: Cavari (2019)
##############################################################################################################################
approval <- read_dta("main-cavari-psq.dta")

# Estimating the Original Model
approval$rr <- ifelse(approval$PRESIDENT == 40, 1, 0)
approval$ghwb <- ifelse(approval$PRESIDENT == 41, 1, 0)
approval$bc <- ifelse(approval$PRESIDENT == 42, 1, 0)
approval$gwb <- ifelse(approval$PRESIDENT == 43, 1, 0)
approval$bo <- ifelse(approval$PRESIDENT == 44, 1, 0)

approval$ECONAPP_ECONMIP <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS
approval$FPAPP_FPMIP <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN

# Create matrix of x regressors
approval.small <- approval %>%
				dplyr::select(APPROVE, APPROVE_ECONOMY, APPROVE_FOREIGN, MIP_MACROECONOMICS, MIP_FOREIGN, ECONAPP_ECONMIP, FPAPP_FPMIP, 
					APPROVE_L1, PARTY_IN, PARTY_OUT, UNRATE, DIVIDEDGOV, ELECTION, HONEYMOON, ghwb, bc, gwb, bo)

cavari.model <- lm(as.formula(paste0("APPROVE ~ ", paste(names(approval.small), collapse = "+"))), data = approval)

summary(cavari.model)
bgtest(cavari.model, order = 5, type = "F")

## Extend to the general model

# Creating Lags 
approval$MIP_FOREIGN_L1 <- lshift(approval$MIP_FOREIGN, 1)
approval$MIP_MACROECONOMICS_L1 <- lshift(approval$MIP_MACROECONOMICS, 1)
approval$APPROVE_ECONOMY_L1 <- lshift(approval$APPROVE_ECONOMY, 1)
approval$APPROVE_FOREIGN_L1 <- lshift(approval$APPROVE_FOREIGN, 1)
approval$UNRATE_L1 <- lshift(approval$UNRATE, 1)
approval$PARTY_IN_L1 <- lshift(approval$PARTY_IN, 1)
approval$PARTY_OUT_L1 <- lshift(approval$PARTY_OUT, 1)

approval$ECONAPP_L1_ECONMIP_L1 <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_L1
approval$FPAPP_L1_FPMIP_L1 <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_L1

approval$ECONAPP_ECONMIP_L1 <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_L1
approval$FPAPP_FPMIP_L1 <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_L1

approval$ECONAPP_L1_ECONMIP <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS
approval$FPAPP_L1_FPMIP <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN

#Testing stationarity of all variables and interactions
df_ECONAPP_ECONMIP <- ur.df(approval$ECONAPP_ECONMIP, type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP)
df_FPAPP_FPMIP <- ur.df(approval$FPAPP_FPMIP, type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_FPMIP)

df_ECONAPP_ECONMIP_L1 <- ur.df(approval$ECONAPP_ECONMIP_L1[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP_L1)
df_FPAPP_FPMIP_L1 <- ur.df(approval$FPAPP_FPMIP_L1[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_FPMIP_L1)

df_ECONAPP_L1_ECONMIP <- ur.df(approval$ECONAPP_L1_ECONMIP[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_ECONAPP_ECONMIP)
df_FPAPP_L1_FPMIP <- ur.df(approval$FPAPP_L1_FPMIP[2:140], type = "drift", selectlags = "AIC")
interp_urdf(df_FPAPP_L1_FPMIP)

df_APPROVAL <- ur.df(approval$APPROVE, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVAL)

df_APPROVE_ECONOMY <- ur.df(approval$APPROVE_ECONOMY, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVE_ECONOMY)
df_APPROVE_FOREIGN <- ur.df(approval$APPROVE_FOREIGN, type = "drift", selectlags = "AIC")
interp_urdf(df_APPROVE_FOREIGN)
df_MIP_MACROECONOMICS <- ur.df(approval$MIP_MACROECONOMICS, type = "drift", selectlags = "AIC")
interp_urdf(df_MIP_MACROECONOMICS) #Unit root!!!
df_MIP_FOREIGN <- ur.df(approval$MIP_FOREIGN, type = "drift", selectlags = "AIC")
interp_urdf(df_MIP_FOREIGN) #Unit root!!!
df_PARTY_IN <- ur.df(approval$PARTY_IN[1:138], type = "drift", selectlags = "AIC")
interp_urdf(df_PARTY_IN) #Unit root!!!
df_PARTY_OUT <- ur.df(approval$PARTY_OUT[1:138], type = "drift", selectlags = "AIC")
interp_urdf(df_PARTY_OUT)
df_UNRATE <- ur.df(approval$UNRATE, type = "drift", selectlags = "AIC")
interp_urdf(df_UNRATE) #Unit root!!!
df_DIVIDEDGOV <- ur.df(approval$DIVIDEDGOV, type = "drift", selectlags = "AIC")
interp_urdf(df_DIVIDEDGOV)

###Difference Unit Roots
approval$UNRATE_D <- approval$UNRATE - approval$UNRATE_L1
approval$MIP_MACROECONOMICS_D <- approval$MIP_MACROECONOMICS - approval$MIP_MACROECONOMICS_L1
approval$MIP_FOREIGN_D <- approval$MIP_FOREIGN - approval$MIP_FOREIGN_L1

# Create lags
approval$MIP_FOREIGN_D_L1 <- lshift(approval$MIP_FOREIGN_D, 1)
approval$MIP_MACROECONOMICS_D_L1 <- lshift(approval$MIP_MACROECONOMICS_D, 1)
approval$UNRATE_D_L1 <- lshift(approval$UNRATE_D, 1)

# Create lagged and cross-lagged interactions for differenced variable
# t/t
approval$ECONAPP_ECONMIP_D <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_D
approval$FPAPP_FPMIP_D <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_D

# t-1/t-1
approval$ECONAPP_L1_ECONMIP_D_L1 <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_D_L1
approval$FPAPP_L1_FPMIP_D_L1 <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_D_L1

# t/t-1
approval$ECONAPP_ECONMIP_D_L1 <- approval$APPROVE_ECONOMY*approval$MIP_MACROECONOMICS_D_L1
approval$FPAPP_FPMIP_D_L1 <- approval$APPROVE_FOREIGN*approval$MIP_FOREIGN_D_L1

# t-1/t
approval$ECONAPP_L1_ECONMIP_D <- approval$APPROVE_ECONOMY_L1*approval$MIP_MACROECONOMICS_D
approval$FPAPP_L1_FPMIP_D <- approval$APPROVE_FOREIGN_L1*approval$MIP_FOREIGN_D

approval_standard <- approval %>%
						filter(!is.na(UNRATE_D_L1), !is.na(PARTY_IN))

#  Report these in the same way as the J_J replication in the Supplemental Materials
adf.test(approval_standard$APPROVE_ECONOMY)$p.value
adf.test(approval_standard$APPROVE_FOREIGN)$p.value
adf.test(approval_standard$MIP_MACROECONOMICS)$p.value
adf.test(approval_standard$MIP_FOREIGN)$p.value
adf.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
adf.test(approval_standard$MIP_FOREIGN_D)$p.value
adf.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
adf.test(approval_standard$FPAPP_FPMIP_D)$p.value
adf.test(approval_standard$PARTY_IN)$p.value
adf.test(approval_standard$PARTY_OUT)$p.value
adf.test(approval_standard$UNRATE)$p.value
adf.test(approval_standard$UNRATE_D)$p.value
adf.test(approval_standard$DIVIDEDGOV)$p.value
adf.test(approval_standard$ELECTION)$p.value
adf.test(approval_standard$HONEYMOON)$p.value


Box.test(approval_standard$APPROVE_ECONOMY)$p.value
Box.test(approval_standard$APPROVE_FOREIGN)$p.value
Box.test(approval_standard$MIP_MACROECONOMICS)$p.value
Box.test(approval_standard$MIP_FOREIGN)$p.value
Box.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
Box.test(approval_standard$MIP_FOREIGN_D)$p.value
Box.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
Box.test(approval_standard$FPAPP_FPMIP_D)$p.value
Box.test(approval_standard$PARTY_IN)$p.value
Box.test(approval_standard$PARTY_OUT)$p.value
Box.test(approval_standard$UNRATE)$p.value
Box.test(approval_standard$UNRATE_D)$p.value
Box.test(approval_standard$DIVIDEDGOV)$p.value
Box.test(approval_standard$ELECTION)$p.value
Box.test(approval_standard$HONEYMOON)$p.value


kpss.test(approval_standard$APPROVE_ECONOMY)$p.value
kpss.test(approval_standard$APPROVE_FOREIGN)$p.value
kpss.test(approval_standard$MIP_MACROECONOMICS)$p.value
kpss.test(approval_standard$MIP_FOREIGN)$p.value
kpss.test(approval_standard$MIP_MACROECONOMICS_D)$p.value
kpss.test(approval_standard$MIP_FOREIGN_D)$p.value
kpss.test(approval_standard$ECONAPP_ECONMIP_D)$p.value
kpss.test(approval_standard$FPAPP_FPMIP_D)$p.value
kpss.test(approval_standard$PARTY_IN)$p.value
kpss.test(approval_standard$PARTY_OUT)$p.value
kpss.test(approval_standard$UNRATE)$p.value
kpss.test(approval_standard$UNRATE_D)$p.value
kpss.test(approval_standard$DIVIDEDGOV)$p.value
kpss.test(approval_standard$ELECTION)$p.value
kpss.test(approval_standard$HONEYMOON)$p.value

cavari.model.differenced <- lm(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + 
                     MIP_FOREIGN_D + ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + 
                     APPROVE_L1 + PARTY_IN + PARTY_OUT + UNRATE_D + 
                     DIVIDEDGOV + ELECTION + HONEYMOON + ghwb + bc + gwb + bo, data = approval_standard)
summary(cavari.model.differenced)

# Following our own advice, we will test to see whether any parameters can be safely omitted via LASSO
polywog.fit <- polywog(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + MIP_FOREIGN_D +
                         ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + 
                         APPROVE_ECONOMY_L1 + APPROVE_FOREIGN_L1 + MIP_MACROECONOMICS_D_L1 + MIP_FOREIGN_D_L1 +
                         ECONAPP_L1_ECONMIP_D_L1 + FPAPP_L1_FPMIP_D_L1 +
                         ECONAPP_ECONMIP_D_L1 + FPAPP_FPMIP_D_L1 + ECONAPP_L1_ECONMIP_D + FPAPP_L1_FPMIP_D +
                         APPROVE_L1 + PARTY_IN + PARTY_OUT + UNRATE_D + 
                         PARTY_IN_L1 + PARTY_OUT_L1 + UNRATE_D_L1 + 
                         DIVIDEDGOV + ELECTION + HONEYMOON +
                         ghwb + bc + gwb + bo,
                       data = approval_standard, unpenalized = c("ghwb", "bc", "gwb", "bo"),
                method = "alasso", family = "gaussian", degree = 1, nfolds = 10)
summary(polywog.fit)

# Interactive structure implied from polywog. A mess of cross-lags
#             APP          MIP_D           INTAXN
# ECON        t, t-1	   no, no			 t/t-1, t-1/t
# F           t, t-1       no, no          t-1/t, t-1/t-1

#Refit with selected parameters. we will keep the general structure but pare down the other X
polywog.model <- lm(APPROVE ~ APPROVE_ECONOMY + APPROVE_FOREIGN + MIP_MACROECONOMICS_D + MIP_FOREIGN_D + # all t terms
                  APPROVE_ECONOMY_L1 + APPROVE_FOREIGN_L1 + MIP_MACROECONOMICS_D_L1 + MIP_FOREIGN_D_L1 + # all t-1 terms
                  ECONAPP_ECONMIP_D + FPAPP_FPMIP_D + # t/t
                  ECONAPP_L1_ECONMIP_D_L1 + FPAPP_L1_FPMIP_D_L1 + #t-1/t-1
                  ECONAPP_ECONMIP_D_L1 + ECONAPP_L1_ECONMIP_D + #t/t-1 and t-1/t 
                  FPAPP_FPMIP_D_L1 + FPAPP_L1_FPMIP_D + #t/t-1 and t-1/t
					APPROVE_L1 + PARTY_IN + PARTY_OUT_L1 +
                    ghwb + bc + gwb + bo, data = approval_standard)

summary(polywog.model)

## Function for returning instantaneous and long-run effects over range of MIP_macro (z)
effects.mipmacro <- function(change.MIP_MACROECONOMICS_D) {
  thetas <- coef(polywog.model)
  varcov <- vcov(polywog.model)
  theta.tilde <- data.frame(mvrnorm(5000,thetas,varcov))

	theta.tilde <- theta.tilde[,which(names(theta.tilde) %in% c("APPROVE_L1", "APPROVE_ECONOMY", "MIP_MACROECONOMICS_D", 
																"APPROVE_ECONOMY_L1", "MIP_MACROECONOMICS_D_L1",
																"ECONAPP_ECONMIP_D", "ECONAPP_L1_ECONMIP_D_L1",
																"ECONAPP_ECONMIP_D_L1", "ECONAPP_L1_ECONMIP_D"))]
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_L1")] <- "alpha_1"
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_ECONOMY")] <- "beta_0"   # x_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_ECONOMY_L1")] <- "beta_1"		# x_{t-1} or lag x
	colnames(theta.tilde)[(colnames(theta.tilde) == "MIP_MACROECONOMICS_D")] <- "beta_2"      # z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "MIP_MACROECONOMICS_D_L1")] <- "beta_3"     # z_{t-1}
	colnames(theta.tilde)[(colnames(theta.tilde) == "ECONAPP_ECONMIP_D")] <- "beta_4"       # x_t z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "ECONAPP_L1_ECONMIP_D")] <- "beta_5"			 # x_{t-1} z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "ECONAPP_ECONMIP_D_L1")] <- "beta_6"				 # x_t z_{t-1}
	colnames(theta.tilde)[(colnames(theta.tilde) == "ECONAPP_L1_ECONMIP_D_L1")] <- "beta_7"				 # x_{t-1} z_{t-1}

	# mip <- seq(unname(summary(approval_standard$MIP_MACROECONOMICS_D)[1]),
    #              unname(summary(approval_standard$MIP_MACROECONOMICS_D)[6]), length.out = 250) # vary level
	# wacky range. let's just do normal from -10 to 10
	mip <- seq(-10, 10, length.out = 250)
	
	# create hypothetical values based on these
	hyp_zt <- mip  # range of level of changes
	hyp_zt_l1 <- mip - change.MIP_MACROECONOMICS_D # level - change applied = lag (level before change)
	hyp_zt_p1 <- mip + change.MIP_MACROECONOMICS_D # level + change applied = future (level after change)

	# create list to loop over
	hyplist <- as.list(NULL)
	# create storage for results
	effects <- data.frame(matrix(ncol=7,nrow=length(mip)))

	colnames(effects) <- c("mip","inst.lo","inst.est.median","inst.hi","total.lo","total.est.median","total.hi")

	for(i in 1:length(mip)) {  # for all of the values of MIP
		### create temp storage
		hyplist[[i]] <- theta.tilde

		## calculate elements from instantaneous effect and long-run effect formulae
		# instantaneous
		# beta_0 already an element of hyplist
		hyplist[[i]]$beta4zt <- hyplist[[i]]$beta_4*hyp_zt[i]   # beta_4 $\times$ the range of hypothetical z_ts
		hyplist[[i]]$beta6ztminus1 <- hyplist[[i]]$beta_6*hyp_zt_l1[i] # beta_6 $\times$ the post-change z_t (z_{t+1})
		
		# LRE
		# beta_0 already an element of hyplist
		# beta_1 already an element of hyplist
		# beta_4 * z_t already an element of hyplist (from instantaneous)
		hyplist[[i]]$beta5ztplus1 <- hyplist[[i]]$beta_5*hyp_zt_p1[i]
		# beta_6 * z_{t-1} already an element of hyplist (from instantaneous)
		hyplist[[i]]$beta7zt <- hyplist[[i]]$beta_7*hyp_zt[i]   # beta_7 $\times$ the range of hypothetical z_ts
		
		### calculate quantities of interest
    	hyplist[[i]]$inst.effect <- rowSums(hyplist[[i]][,c("beta_0", "beta4zt","beta6ztminus1")])
    	hyplist[[i]]$total.effect <- (rowSums(hyplist[[i]][,c("beta_0", "beta_1","beta4zt","beta5ztplus1", "beta6ztminus1","beta7zt")]))/ # num
                                                          					(1 - hyplist[[i]]$alpha_1)								# denom
		effects$mip[i] <- hyp_zt[i]
		effects$inst.lo[i] <- unname(quantile(hyplist[[i]]$inst.effect, .025))
		effects$inst.est.median[i] <- unname(summary(hyplist[[i]]$inst.effect)["Median"])
		effects$inst.hi[i] <- unname(quantile(hyplist[[i]]$inst.effect, .975))
		effects$total.lo[i] <- unname(quantile(hyplist[[i]]$total.effect, .025))
		effects$total.est.median[i] <- unname(summary(hyplist[[i]]$total.effect)["Median"])
		effects$total.hi[i] <- unname(quantile(hyplist[[i]]$total.effect, .975))
	}
	return(effects)
}

effects.econ <- effects.mipmacro(change.MIP_MACROECONOMICS_D = mean(approval_standard$MIP_MACROECONOMICS_D))

# plot
cavari.econ.inst <- ggplot(effects.econ, aes(x = mip, y = inst.est.median)) + 
  theme_bw() + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  geom_ribbon(aes(ymin = inst.lo, ymax = inst.hi), alpha=.25) +
  geom_line() +
  geom_hline(aes(yintercept = 0), color = "gray30", linetype = "dashed") +
  # coord_cartesian(ylim=c(-7, 7)) +
  labs(title = "Instantaneous effects\n(economic approval)", x = expression(paste(Delta, "Salience of economy as MIP")), 
  	y = expression(paste(Delta, "Presidential approval"))) # +
 # theme(text=element_text(family="minpro"))
 

cavari.econ.total <- ggplot(effects.econ, aes(x = mip, y = total.est.median)) + 
  theme_bw() + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  geom_ribbon(aes(ymin = total.lo, ymax = total.hi), alpha=.25) +
  geom_line() +
  geom_hline(aes(yintercept = 0), color = "gray30", linetype = "dashed") +
  # coord_cartesian(ylim=c(-7, 7)) +
  labs(title = "Total effects\n(economic approval)", x = expression(paste(Delta, "Salience of economy as MIP")), 
  	y = expression(paste(Delta, "Presidential approval"))) # +
 # theme(text=element_text(family="minpro"))




## Function for returning instantaneous and long-run effects over range of MIP_fp (z)
effects.mipfp <- function(change.MIP_FOREIGN_D){
	thetas <- coef(polywog.model)
	varcov <- vcov(polywog.model)
	theta.tilde <- data.frame(mvrnorm(5000,thetas,varcov))

	theta.tilde <- theta.tilde[,which(names(theta.tilde) %in% c("APPROVE_L1", "APPROVE_FOREIGN", "MIP_FOREIGN_D", 
																"APPROVE_FOREIGN_L1", "MIP_FOREIGN_D_L1",
																"FPAPP_FPMIP_D", "FPAPP_L1_FPMIP_D_L1",
																"FPAPP_FPMIP_D_L1", "FPAPP_L1_FPMIP_D"))]
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_L1")] <- "alpha_1"
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_FOREIGN")] <- "beta_0"   # x_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "APPROVE_FOREIGN_L1")] <- "beta_1"		# x_{t-1} or lag x
	colnames(theta.tilde)[(colnames(theta.tilde) == "MIP_FOREIGN_D")] <- "beta_2"      # z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "MIP_FOREIGN_D_L1")] <- "beta_3"     # z_{t-1}
	colnames(theta.tilde)[(colnames(theta.tilde) == "FPAPP_FPMIP_D")] <- "beta_4"       # x_t z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "FPAPP_L1_FPMIP_D")] <- "beta_5"			 # x_{t-1} z_t
	colnames(theta.tilde)[(colnames(theta.tilde) == "FPAPP_FPMIP_D_L1")] <- "beta_6"				 # x_t z_{t-1}
	colnames(theta.tilde)[(colnames(theta.tilde) == "FPAPP_L1_FPMIP_D_L1")] <- "beta_7"				 # x_{t-1} z_{t-1}

	# mip <- seq(unname(summary(approval_standard$MIP_FOREIGN_D)[1]),
    #              unname(summary(approval_standard$MIP_FOREIGN_D)[6]), length.out = 250) # vary level
	# wacky range. let's just do normal from -10 to 10
	mip <- seq(-10, 10, length.out = 250)
	
	
	# create hypothetical values based on these
	hyp_zt <- mip  # range of level of changes
	hyp_zt_l1 <- mip - change.MIP_FOREIGN_D # level - change applied = lag (level before change)
	hyp_zt_p1 <- mip + change.MIP_FOREIGN_D # level + change applied = future (level after change)

	# create list to loop over
	hyplist <- as.list(NULL)
	# create storage for results
	effects <- data.frame(matrix(ncol = 7, nrow = length(mip)))

	colnames(effects) <- c("mip", "inst.lo", "inst.est.median", "inst.hi", "total.lo", "total.est.median", "total.hi")

	for(i in 1:length(mip)) {  # for all of the values of MIP
		### create temp storage
		hyplist[[i]] <- theta.tilde

		## calculate elements from instantaneous effect and long-run effect formulae
		# instantaneous
		# beta_0 already an element of hyplist
		hyplist[[i]]$beta4zt <- hyplist[[i]]$beta_4*hyp_zt[i]   # beta_4 $\times$ the range of hypothetical z_ts
		hyplist[[i]]$beta6ztminus1 <- hyplist[[i]]$beta_6*hyp_zt_l1[i] # beta_6 $\times$ the post-change z_t (z_{t+1})
		
		# LRE
		# beta_0 already an element of hyplist
		# beta_1 already an element of hyplist
		# beta_4 * z_t already an element of hyplist (from instantaneous)
		hyplist[[i]]$beta5ztplus1 <- hyplist[[i]]$beta_5*hyp_zt_p1[i]
		# beta_6 * z_{t-1} already an element of hyplist (from instantaneous)
		hyplist[[i]]$beta7zt <- hyplist[[i]]$beta_7*hyp_zt[i]   # beta_7 $\times$ the range of hypothetical z_ts
		
		### calculate quantities of interest
    	hyplist[[i]]$inst.effect <- rowSums(hyplist[[i]][,c("beta_0", "beta4zt","beta6ztminus1")])
    	hyplist[[i]]$total.effect <- (rowSums(hyplist[[i]][,c("beta_0", "beta_1","beta4zt","beta5ztplus1", "beta6ztminus1","beta7zt")]))/ # num
                                                          					(1 - hyplist[[i]]$alpha_1)								# denom
		effects$mip[i] <- hyp_zt[i]
		effects$inst.lo[i] <- unname(quantile(hyplist[[i]]$inst.effect, .025))
		effects$inst.est.median[i] <- unname(summary(hyplist[[i]]$inst.effect)["Median"])
		effects$inst.hi[i] <- unname(quantile(hyplist[[i]]$inst.effect, .975))
		effects$total.lo[i] <- unname(quantile(hyplist[[i]]$total.effect, .025))
		effects$total.est.median[i] <- unname(summary(hyplist[[i]]$total.effect)["Median"])
		effects$total.hi[i] <- unname(quantile(hyplist[[i]]$total.effect, .975))
	}
	return(effects)
}

effects.fp <- effects.mipfp(change.MIP_FOREIGN_D = mean(approval_standard$MIP_FOREIGN_D))

# plot
cavari.fp.inst <- ggplot(effects.fp, aes(x = mip, y = inst.est.median)) + 
  theme_bw() + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  geom_ribbon(aes(ymin = inst.lo, ymax = inst.hi), alpha=.25) +
  geom_line() +
  geom_hline(aes(yintercept = 0), color = "gray30", linetype = "dashed") +
  # coord_cartesian(ylim=c(-7, 7)) +
  labs(title = "Instantaneous effects\n(foreign policy approval)", x = expression(paste(Delta, "Salience of foreign policy as MIP")), 
  	y = expression(paste(Delta, "Presidential approval"))) # +
 # theme(text=element_text(family="minpro"))


cavari.fp.total <- ggplot(effects.fp, aes(x = mip, y = total.est.median)) + 
  theme_bw() + 
  theme(panel.grid.minor = element_blank(), panel.border = element_rect(fill = NA, size = 1, linetype = "solid")) +
  geom_ribbon(aes(ymin = total.lo, ymax = total.hi), alpha=.25) +
  geom_line() +
  geom_hline(aes(yintercept = 0), color = "gray30", linetype = "dashed") +
  # coord_cartesian(ylim=c(-7, 7)) +
  labs(title = "Total effects\n(foreign policy approval)", x = expression(paste(Delta, "Salience of foreign policy as MIP")), 
  	y = expression(paste(Delta, "Presidential approval"))) # +
 # theme(text=element_text(family="minpro"))



######### Figure 5 #########
pdf("fig5.pdf", width = 6, height = 6)
multiplot(cavari.econ.inst, cavari.fp.inst, cavari.econ.total, cavari.fp.total, cols = 2)
dev.off()


